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NUMERICAL  SOLUTION  OF  A  MODEL  PROBLEM  FROM  COLLAPSE  LOAD  ANALYSIS 


Michael  L.  Overton 

Courant  Institute  of  Mathematical  Sciences 
New  York  University 
New  York,  N.  Y.  10012 
U.S.A. 


We  consider  a  model  problem  from  collapse  load  analysis 
discussed  recently  by  Strang.  The  analytic  solution  is 
a  characteristic  function  with  a  jump  discontinuity.  We 
develop  a  method  for  solving  a  discretized  version  of  the 
model  problem,  which  requires  the  minimization  of  a  convex 
piecewise  differentiable  function.  Numerical  results  are 
presented. 


INTRODUCTION 


In  collapse  load  analysis  the  following  problem  arises  (see  Strang  [16]) 


min    IvuB 
u  -'n 

subject  to  the  Dirichlet  boundary  conditions 

u(x,y)  =  f(x,y)  on  an 


(l.Ta) 


(1.1b) 


and  the  constraint 


c(x,y)  u(x,y)  =  1 


(1.1c) 


Q 


Here  llvull  denotes  the  E 
We  shall  restrict  Q.  to 
usually  smooth  function 
straint  (1.1c)  does  not 
solution  if,  for  exampl 
appropriate  space  of  fu 
tion.  Strang  shows,  us 
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inside  r  and  zero  value 
curve  and  the  broken  li 
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e,  f(x,y)  E  0.  In 
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ing  a  result  of  Fl 
n  u(x,y)  is  a  chara 
n  Q.  Specifically, 
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nes  represent  an: 


e  gradient  of  u(x,y),  i.e. 
de  tv/o.  The  functions   f   and-'  c  are 
d  n  respectively.  Sometimes  the  con- 
s  required  to  obtain  a  nontrivial 
order  for  a  minimum  to  be  achieved,  the 

is  BV,  the  functions  of  bounded  varia- 
eming  [9],  that  when  c(x,y)  =  1, 
cteristic  function  with  a  jump  disconti- 

u(x,y)  has  constant  positive  value 
following  figure,  where  r  is  the  solid 


(1.2) 


j_.  _  ;=~ 


<r J 


The  curve  r  consists  of  four  circular  arcs  of  radius  p  =  1/(1+  ^11)   -   0.530, 
plus  parts  of  the  sides  of  the  square.  Once  it  is  recognized  that  the  solution 
is  a  characteristic  function,  (1.1)  reduces  to  a  classical  isoperimetric  problem, 
since  the  quantity  to  be  minimized  reduces  to  measuring  the  length  of  r  while  the 
constraint  fixes  the  area  inside  r.  The  solution  therefore  defines  the  region 
in  the  square  with  minimal  perimeter  given  fixed  area  (maximal  area  given  fixed 
perimeter) . 

In  this  paper  we  are  concerned  with  Obtaining  the  numerical  solution  of  a  discrete 
approximation  to  (1.1).  Let  us  triangulate  n  as  follows: 


^vr 


0,N-1 


0,0 


.N-1,N-1 


t  ''-J^^l^.. 


Let  h  =  1/(N-1)     be  the  mesh  rire, jwh?r^§,ltf]i^re^§re  N"mesh  points  in  each  direc- 
tion.    We  replace  u  in  (1.1)  by  a  piecewise^li''rie4r  finTte  element  approximation  v, 
obtaining  a  finite  dimensional  optimization  problem.     The  variables  are  the 
unknown  function  values  v(x-j  ,y-;) ,  at  the-giyen  mesh  points,   i   =  0,...,N-1, 
j  =  0,...,N-1.     Because  v  is  pTecewiserli near, 'the  norm  of  its  gradient 


=  (v  ,v  )''"  on  a  "lower"  triangle 
X  y 


IX. 


-I'^j^ 


(Xi.ryj-i) 


I  > 


(Xi.yj.i) 


multiplied  by  the  area  of  the  triangle,  is  given  by 


Br.  .J  = 
ijl 


Let  us  write 


r-h 


2  (v(Xi,yj.T)-v(x._^,yj_T)) 
L|(v(x,..pyj)-v(x..pyj..i))J 


r.jT  =  A.  .^  V  E  IR^ 

N  2 

where  v  is  the  vector  e  IR   of  unknowns.  The  matrix  A- .,  has  dimension  N  ^  2 


lA         "tiT      "ti-^     uimciijiuii     \\ 

and  has  only  two  nonzero  components  in  each  column.  We  use  r- -^  and  A.jj2  f 
the  corresponding  notation  for  "upper"  triangles.  Thus  the  discretized  pro 
is:  i 


or 
problem 


min 


r ...  II 
ijk 


V     i,j,k 
(i=l,...,N-l,  j=l,...,N-l,  k=l,2) 


(1.3a) 


-3- 
subject  to  the  boundary  conditions 

v(x.  ,yj)  =  ■f(x.,yj)  for  (x-.y^)  e  sn  (1.3b) 

and  the  constraint  •     " 

-V    I     c(x..y.),y(x.,y.)  =  1    .  (1.3c) 

N"^  i,j     .  ^  ,:--^    ^ 

Although   (1.3)   is  'i  convex  finite-dimensional   optimization  problem, 
it  is  not  easily  solved,  because  the  objective  function  to  be  minimized  is  only 
piecewise  differentiable  with-respefct~~tl5"y^'rTmy  of  the  r-jii^  equal   zero.   Now  the 
residual     r-jj^^  is  zero  precisely  if  the  sCluti-on  is  constant  on  the  corresponding 
triangle.     The  solution  of  (1.1)    (with  c-r—l- -,•  f •  =  0)     is  constant  everywhere 
except  across    r    where  it  jumps.     We  expect  ;tl>eref ore  that  the  solution  to  the 
discretized  problem  will   be  constartt^v^r.  iiij;isjLM:^.  but  varying  in  some  band 
around     r.     In      this  case  (1.3)  is;  indfeed  a  nohdifferentiable  optimization 
problem.  'J^    '  ''' 


./    / 

AN  ALGORITHM  FOR  MINIMIZING  'A-^SUM  CfTElfcn 

We  have  recently  proposed.^4^  SQ ^a]^goH thiji ^f or  solving  nondifferentiable  optimi- 
zation problems  of  this  fdrnj...  Soi^WhSt 'ratf re  ^eneratly,  it  solves  the  convex 

.:/;::intn::F(v)^^:f^fr[(y|r-  (2.1) 

.  vem"  ^-^   -.,1   ...   - 

where  - 

r^-(v)  =  aJv  -  b^.  €   :IR^. 

We  are  changing  the  notation  slightly,  writing  i  where  before  we  wrote  the  triple 
(i,j,k),  and  now  allowing  the  residual  r-j  to  be  affine  rather  than  linear  in  v. 
As  before,  II -i  denotes  the  Euclidean  vector  norm.  Linear  constraints  may  be 
included  without  difficulty,  but  we  omit  them  for  convenience.  The  problem  (2.1) 
dates  back  to  Fermat  and  has  an  interesting  history.  In  its  simplest  form,  with 
n  =  £.  =  2,  it  asks  for  the  point  in  the  plane  which  minimizes  the  sum  of 
distances  betv;een  it  and  m  given  points.  A  more  interesting  version  involves  a 
weighted  sum  of  the  distances.  With  this  variation,  the  solution  coincides  with 
one  of  the  given  points  if  the  corresponding  weight  is  large  enough.  This  problem 
is  also  associated  with  the  names  of  Steiner  and  Weber  and  is  often  known  as  the 
single  facility  location  problem.  The  multi facility  location  problem  arises 
when  more  than  one  point  in  the  plane  is  to  be  determined,  i.e.  n  >  i  =  2.  The 
weighted  sum  of  norms  to  be  minimized  then  involves  the  distances  between  each 
pair  of  unknowns  as  well  as  the  distances  between  the  unknown  and  fixed  points. 
This  particular  problem  is  discussed  in  more  detail  by  Calamai  and  Conn  [2,3,4], 
who  give  an  algorithm  closely  related  to  ours. 

The  algorithm  as  described  in  [14]  is  intended  to  solve  (2.1)  when  the  matrices 
{Ai}  are  dense  ana  n  is  not  too  large;  a  full  rank  condition  is  also  required. 
None  of  these  conaitions  holds  in  (1.3).  The  purpose  of  this  paper  is  to  explain 
how  to  adapt  the  algorithm  to  solve  (1.3),  and  to  report  on  the  advantages  and 
difficulties  of  such  an  approach. 

Notation.  Using  ":  to  denote  gradient  with  respect  to  v,  we  have: 

9i(^)  =  '"^(^)"  =  177(771  ^^■(^^ 
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1  ^  '        I  r .  ( V )  0 

provided  that  llr.  (v)l!  ^  0.  Note  that  the  Hessian  term  G^-(v)  is  unbounded  as 
Dr-(v)ll  -»-  0,  but  that  the  gradient  term  g-j{v)  remains  bounded,  although  it  is  of 
course  discontinuous  at  v  if  r^{v)  =  0.   Now  let  us  define 

g(v)  *    I     -^  g.(v) 

and  •,  ''i^m::,^.  ' 

G(v)  =    I  G.(v)  . 

lr.{v)nj«0  ' 

If  the  gradient  and  Hessian  of  F(v)  are  defined  it  is  clear  that  they  are  given  by 
g(v)  and  G(v)   respectively.     The  discontinuities  in  the  gradient  and  consequent 
unboundedness  in  the  Hessian  are  what  cause  the  difficulty  in  solving   (2.1). 

In  order  to  be  able  to  handle  the  dfscentihuiti'es- consider  the  following  idea. 
Define  the  active  set  at  a  point  v'as-'-'     ''     "'■"-  s"*'    ^   -• 


J{vO-^in'^  trt(v9rs'0}V  .v^  (2.2) 

i.e.   those  indices  associated  w'ith-t&rb  r%sidtial4'.- -'We':  refer  to     r-  as  an  active 
residual   or  active  term  at  v  if  i~e'~Cl(vO-.-   Thfe-'1<jea  Of:the  algorithm  is  to  project 
the  objective  function  F(v)   into-  the^ ^'1  near' frtah4 fold' 5/here  zero  residuals  remain 
unchanged.     In  this  space  we  see- that"  F(v)"  Is  lt>cal-1y  tbntinuously  differentiable, 
since  discontinuities  occur  only,  in  di recti orvs  cy'ossijig  the  manifold.     To  make 
this  precise  let  "    '- -'f  --' •'  ■"  f  "    "r.=  ?r 

^  ^  '        ■      '       ' 

■  A  =  A(v)  =  [A.  A.  ...]  where  J  =  {i, .i,,.. .} 
M  ^'2-''    "     -      ■  '  ^ 

i.e.  the  matrix  whose  columns  are  the  "^coe'ffi-cient  matrices  of  the  active  residuals. 
Let  Z  =  Z(v)  be  a  matrix  with  maximal  coTumn  rank  such  that  a'''z  =  0,  i.e.  such 
that  the  columns  of  Z  span  the  null  space  of  A'.  The  significance  of  the  matrix 
Z  is  that 

r^(v+p)  =  r^(v)  =  0  if  i  e   J(v)  and  p  g  R(Z)  , 

the  range  space  of  Z.  Now  for  any  v  consider  the  matrix  Z  and  define  F  restricted 
to  the  space  v  9  R(Z)  as 

F2(P2)  =  F(v  +  Zp^). 

Consider  the  gradient  and  Hessian  of  F^,  which  will  be  called  respectively  the 
projected  gradient  and  projected  Hessian.  Because  the  active  terms  are  fixed  for 
all  p^  ,  we  have 

vF2(P2)  =  Z^g(v+Zp2) 

v^^^iPj)  =   Z^G(v+Zp2)Z  . 

Thus  the  projected  gradient  and  projected  Hessian  are  defined  and  differentiable 
in  a  neighborhood  of  v  regardless  of  whether  any  residuals  are  zero. 

Optimal  ity  conditi  c^.s .  Let  us  consider  conditions  for  a  ooint  v  to  oe  a  solution 
to  (2.1 ), "making  uSe  of  Lhe  above  notation.  Clearly  a  necessary  condition  for  v 
to  be  a  solution  is  that  the  projected  gradient  is  zero,  since  otherwise  F  could 
be  decreased  along  a  direction  in  R(Z(v)).  Thus  for  optimal ity  we  require 
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z"^g(v)  =  0  .  (2.3) 

By  the  definition  of  Z  an  equivalent  condition  is 

g(x)  =  ,  I  Lt.  (2.4) 

l€J(v)  :'  ' 

for  some  vectors  t^-  s  R   for  each;;i  s  J(v).  We  call  the  {t^-}  Lagrange  vectors. 
By  considering  the  change  in  F  along'  d'i-i'ections  not  in  R(Z),  it  can  be  shown 
that  a  necessary  and  sufficient  condition  for  v  to  solve  (2.1)  is  that  (2.4)  hold 
,for  some  vectors  {t^-}  satisfying   -:    ; 

It.l  <  1  for  an  i  e  J(v)  (2.5) 

(see  [12],  [14),  n9]).        .  ?  -  -^v'^-- 

The  algorithm.     The  basic  idea  of  the  algorithm  is  as  follows.     Given  a  point  v^ 
with      an       active     set     J.{vi^);),- we -proceed  to  minimize  the  function  restricted 
to  v^*^^  ®     R(Z(v'^^)),     where  the  zero  residuals  remain  unchanged.     At  the  kth 
step  of  the  iteration,  a  direction  of  search  is  computed  by  solving  the  projected 
Newton  system  of  equations,  utilizing-,    the  projected  gradient  Z^g     and  projected 
Hessian  Z'GZ.     A  special   line  search  is  use'd  to  obtain  the  next  iterate  v(k+1). 
During  the  course  of  thisjit^ratien^wiei-may  r-ediice  other  residuals  to  zero,   in 
which  case  they  are  added- to"  the^, active  set  and^the  restricting  manifold  is 
reduced  in  dimension.     Once- , the' ;^roj^t€d  gradaeflt- z'^g  is  sufficiently  small,  the 
Lagrange  vectors  {t-j}  should  be  £spputed.     Theie:; are  uniquely  defined  if  A  has 
full   column     rank.     If  H  tJI.  <_  il;f-or  a^^l;  ^^e  Or('v).  the  procedure  terminates.    If 
Dtjll>l     for  some  j,  the  j  is  deleted  from  "the  active  set  J,  the  dimension  of  the 
manifold  is  increased,  and  the  process  is  continued. 


C  -1  ■ 


Thedetails  of  the  algorithm  may  be  found  in  [14-3.  ,Here  we  outline  the  main  steps 
which^must  be  performed  to  obtain  a  •new-'-i.-terate  v^'^'*'''  from  v^  '.  It  is  assumed 
that  A  has  full  column  rank  and  that  matrix  factorizations  are  practical;  such 
is  not  the  case  for  (1.3). 

Choose  active  set,  compute  Z,  and  make  almost  active  terms  exactly  active 

if  necessary.  (2.6) 

-  fk) 

Find  the  active  set  J  and  the  matrix  A  at  v^  .  Compute  Z  from  the  QR 

factorization  of  A.  If  some  residuals  are  almost  active,  compute  a  point 

v  which  makes  these  terms  exactly  active.  If  F(v)  <  F(v^'^'),  then 

replace  v^^  by  v  and  restart  this  step;  otherwise  reject  these  terms 

as  inactive. 

Check  optimal ity  (2.7) 

If  the  projected  gradient  norm  HZ  gll^is  small,  compute  the  Lagrange  vectors 
{t.}   using  the  QR  factorization  of  A.  If  Ht^^  <_  Vfor  all  i  s  J  then  stop; 
solution  has  been  approximated.  Otherwise  delete  an  appropriate  term  from  J 
and  set  p  to  a  projected  steepest  descent  direction. 

Solve  projected  Newton  system  (2.8) 

If  IIZ  gll  is  not  small,  solve 

'  (Z^GZ)p,  =  -Z^g 

using  a  Choleski  factorization  and  set  the  direction  of  search  p=  Zp,. 

Line  search  (2.9) 

Use  a  special  line  search,  which  takes  account  of  the  fom  of  F,  to  obtain 
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a  point  v^*^^^^  =  v^^^  +  ap, 


with 


F{v(k'). 


One  interesting  aspect  of  the  algorithm  is  illustrated  by  considering  the  weighted 
Fermat  problem: 


ve  IR 


mi^  Dv-  [  'q  ]e  +  oiOv  -  [  °  ]  D  +  Hv  -  [  J  ] 


(2.10) 


Let  the  initial   iterate  be  given  by,  |ay,  v'^  =    [3,2]    .  We^have  J(vq)  =  0  and 

Z(vq)  =   I.      If     uj  =  1 ,  the  solution^r" — is a  "poi-nt -with  J(v*)  =  0,  and  the 

iterates  generated  by  the  algorithm  cpQverge  quadratically  to  v*  because  of 
standard  properties  of  Newton's  method.     Now  suppose  that     w  =  2.     It  is  easily 
verified     that  the  solution  v*  is   f^rir  >  with  J  containing, one  index.   Now  there 
is  no  obvious  reason  to  suspect  that  the  convergence  of  v^^^   to  v*  will   be  rapid, 
since  F  is  not  differentiable  at  v*  and  the  usual   quadratic  convergence  property 
does  not  apply.      In  fact  the  Hessian  v^F  is  unbounded  in  any  neighborhood  of  v* 
and  undefined  at  v*.      It  is  the  case.,  however,   that  this  ill-conditioning  of  v2f, 
combined  with  the  effect  o|  the  special   Tine     seareh-algorithm,  actually  causes 
quadratic  convergence  to  v   .     Once  'tRe^eonvergeftce'^H  recognized,   the  exact  step 
to  V     may  be  taken  by  (2.6).     Sirice-the-'new'~rest^iet5#ig-aianifold  has  zero  dimen- 
sion, the  Lagrange  vector  is  then: t^mputed- by  (Z-.l^'t  ^ilid-ttie  algorithm  termi- 
nates.    This  quadratic  convergence  reSuH-is- similar^  in  flavor  to  the  well    known 
cubic  convergence  of  the  Rayleigh  quotient  iteration  to  find  an  eigenvector  of  a 
symmetric  matrix.      In  both  cases  the  superlinear  convergence  is  caused  by  the 
fact  that  the  ill-conditioning  of?  the' ma  triced  Vs.  in^  precisely  the  desired  direc- 
tion.    For  further  details   see  •'[14].  -      '-•'-    '  ''r;;-^   :. 


ADAPTING  THE  ALGORITHM  TO  SOLVE"  T^iE  WDD€L- PROBLEM  ''--    'i 

Let  us  now  consider  how  to  adapt  the  algorithm"  of  Section  2  to  solve   (1.3).    It  is 
convenient  to  continue  using  the  notation -oT  Section  2,   i.e.  writing  A-,-   for  A-jj|^, 


n  for  N^,  etc.     The  first  question  is  how 'to  represent  the  active  set  and 
matrix  A.     First  note  that  A  is  certain  to  be  rank     deficient.   This  is  eas 
by  considering  the  following  example: 


1 


the 

ily  seen 


(3.1) 


If  the  solution  v  is  constant  across  triangles  1  and  3,  then  it  must  also  be  con- 
stant across  triangle  2.  Thus  the  three  residuals  being.,zero  are  not  indepen- 
dent properties.  In  general,  the  independent  columns  of  A  can  be  written  in 
the  form 


1 

1 

1 

1 

-1 

-1 

-1 

-1 

1    1 

-1 

■   -1  ■ 

- 

J 

1    1 

,  •:■    rorw  :  .-■■c- 

'  -  J 

• 

^^.. 

(3.2) 


Recall  that  A  p  =  0  is  the  equatign  which  defines  directions  a 
active  set  does  not  change  (zero  re§-iduals  remain  zero).  Thus 
block  of  (3.2)  corresponds  to.-a  regiorv.sucb'as  {3.1 ) ,  where  fi 
v(xT,yj)  are  constrained. to  retain  the- saipg^Gommon  value.  Thi 
however,  free  to  change.-  The- pther^ blocks  01^1(3.2)  correspond 
in  p.  with  constant  v.   r  m-  .-c --;•«'•  -  *■!-''.  o 

Now  let  us  add  columns  to  (3.;2,)  rto  accoun.t  for  the  constraints 
tion  of  search  p  must  also  retain  feasibility.^  §ome  of  the  reg 
solution  may  include  one  or  more  boundary  points,  and  hence  th 
value  must  be  fixed.  Thus  one  column  should  be  added  to  the  c 
making  it  square.  Other  blocks  corresponding  -to  regions  away 
remain  unchanged.  Boundary  points  not  in"  any  of  the  constant 
own  1  X  1  blocks.  Let  us  re-order  the  blocks  so  that  those  cor 
regions  including  boundary  points  appear  first.  The  last  rows 
diagonal  blocks,  correspond  to  free  interior  points.  Finally, 
ponding  to  the  constraint  (1.3c)  must  be  added.  The  matrix  th 


long  which  the 
the  first  diagonal 
ve  function  values 
s  common  value,  is, 
to  other  regions 


,  since  any  direc- 
ions  with  constant 
e  common  function 
orresponding  block, 
from  the  boundary 
regions  form  their 
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Storing  and  factoring  A  would  be  prohibitively  expensive.  Fortunately,  this  is 
not  necessary.  A  full  rank  matrix  Z  which  spans  N(A')  can  be  written  in  the 
following  form: 


Z  = 


(row  a, ) 


(row  Qp) 


(row  y) 


where 


(3.3) 


Note  that  Z  is  not  orthogonal,  Whteh  "WouW  be  desirable  for  numerical  stability. 
However,  this  is'  not  critical.  We  must  rhowever:  ensure  that  we  do  not  divide  by 
a  small  number  in  (3.3).  To  avoid  this';  the'Targest  of  the  values  |c  |,...,|c^| 
is  used,  with  the  rows  interchanged  accordingly.  (Strictly  speaking,  there  might 
be  no  single  free  interior  points,  i.e.  y  >  fi.  but  in  this  unlikely  case  the 
definition  of  Z  can  be  changed  to  have  two  or  more  dense  rows  instead  of  one.) 

The  best  way  to  represent  the  active  set  and  the  matrices  A  and  Z  seems  to  be  as 
follows,  using  a  linked  list  structure.  One  linked  list  is  maintained  for  every 
block  in  A,  including  the  null  blocks  corresponding  to  free  interior  points.  These 
lists  are  of  two  kinds:  fixed  lists,  corresponding  to  square  blocks  in  A,  and 
variable  lists  corresponding  to  other  blocks.  Thus  each  fixed  list  contains 
points  whose  common  function  value  is  fixed  by  the  boundary  conditions,  and  each 
variable  list  contains  points  whose  common  function  value  may  vary.  Initially, 
if  v^'^^  has  no  constant  regions,  each  boundary  point  is  put  in  its  own  fixed  list 
and  each  interior  point  is  put  in  its  own  variable  list.  (Neumann  conditions 
could  also  be  handled  by  putting  two  adjacent  points  from  each  segment  orthogonal 
to  the  boundary  in  their  own  variable  list.)  We  also  maintain  a  table  which 
specifies  the  list  containing  each  point. 

We  now  focus  on  the  best  ways  to  carry  out  steps  (2.6)  through  (2.9)  of  the  algo- 
rithm, using  the  data  structure  just  described.  Let  us  first  consider  (2.6). 
At  the  beginning  of  each  iteration,  the  list  structure  has  ?-fom  which  defines 
the  values  of  J  (and  Z)  at  the  previous  iteration.   It  is  now  first  necessary  to 
add  to  J  any  indices  corresponding  to  residuals  which  became  "exactly"  zero 
during  the  line  searcn.  This  is  done  by  merging  the  appropriate  lists.  For 
example,  consider 
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a 
r 

b 

d 

X 

e 

X 

f 

(3.4) 


supposing  that  at  the  previous  iterati-04v-tl:»ere  was  a  zero  residual  across  tri- 
angle 1  only,  and  that  the  values  of  v  at  points  b,  d,  f  are  fixed  by  the 
Dirichlet  boundary  conditions.  Reflecting  this,  there  would  be  one  fixed  list 
specifying  the  values  of  v  at  (a,b,d),  a  second  fixed  list  specifying  the  value 
of  V  at  f,  and  two  variable  lists  reflecting  the  free  values  at  c  and  e.  Now 
suppose  that  in  the  subsequent  line  search,  residuals  were  reduced  to  zero  across 
triangles  2  and  3  (this  assumes  that  the  ■Oi.r4cWet  data  at  d  and  f  are  the  same, 
and  involves  a  change  only  to  the  value  at  c)'.  The  two  fixed  lists  and  the  vari- 
able list  for  c  would  all  be  merged,  creating  one  fixed  list  and  reflecting  the 
new  definitions  of  J  and  Z.  The  merge  operation  is  trivial  as  it  simply  requires 
appending  one  linked  list  to  another.  —  — 

Step  (2.6)  must  also  consider  forcing  residuals  which  are  "almost"  zero  to  become 
"exactly"  zero.  This  is  a-mor-e-e-offipJ-i-CA-t^d-Ofe-rat-ion  for  two  reasons.  First, 
the  current  list  structure  must  be  saved,  since  the  new  vector  v,  for  which  the 
relevant  residuals  are  exactly  zero,  may  have  a  higher  value  of  F  than  the  current 
iterate,  in  which  case  it  must  be  discarded.  Second,  the  vector  v  must  be  chosen 
to  satisfy  the  constraint  (1.3c).  When  two  variable  lists  are  merged,  with  two 
slightly  different  common  values  of  v,  the  new  variable  common  value  can  be  chosen 
to  ensure  that  (1.3c)  is  satisfied,  provided  that  c(x^-,y-)  is  a  strictly  positive 
function  on  the  whole  mesh.  Thi§  is  the  .case  "for  Strang^s  model  problem  where 
c  =  1.  When  a  variable  and  a  fixed  list  are  merged,  the  variable  common  value 
must  be  changed  to  the  fixed  common  value,  ■:potenti ally  violating  (1.3c).  The 
simplest  way  to  deal  with  this  situation'iSr:to  perform  all  the  merge  operations, 
and  then  restore  feasibility  by  sealing  all  the  variable  values  by  a  constant 
factor.  Again,  it  is  easy  to  see  that  this  is  possible  provided  c  is  a  positive 
function.  The  merging  of  two  fixed  lists  is  not  permitted  when  the  residual  is 
not  "exactly"  zero,  i.e.  within  machine  accuracy. 

It  is  of  some  interest  to  note  that  in  practice,  many  residuals  are  reduced 
exactly  to  zero  by  the  line  search.  This  contrasts  with  the  situation  for  the 
single  and  multi-facility  location  problems,  e.g.  (2.10)  with  oj  =  2,  where  the 
convergence  of  the  residual  to  zero  is  an  iterative  process  with  quadratic  conver- 
gence. The  reason  for  the  different  behavior  here  is  not  the  rank  deficiency  of  A 
per  se,  since  this  is  also  present  in  a  different  form  for  the  multifacility 
location  problem.  Consider  again  the  example  (3.4).  Given  any  p  e  R(Z),  the  only 
component  relevant  to  the  proposition  of  incorporating  triangle  2  into  the 
constant  region  of  triangle  1  is  the  change  in  the  value  of  v  at  c,  or,  if  (a,b,d) 
are  in  a  variable  list,  the  changes  of  the  value  at  c  and  the  value  at  (a,b,d). 
Thus  it  is  normally  possible  to  pick  a  steplength  a  which  will  make  v  +  ap 
constant  across  triangles  1  and  2  (when  (a,b,d)  are  in  a  variable  list,  it  is 
just  a  matter  of  weighting  the  two  values  correctly).  In  other  words,  for  almost 
any  p  e  R(Z),  the  extra  residual  can  be  reduced  to  zero  along  p.  In  the  weighted 
Fermat  problem  (2.10),  with  u  =  2,  almost  no  direction  points  from  a  given  v  in 
the  plane  to  the  solution  v*,  the  only  point  where  the  relevant  residual  is  zero. 

Let  us  now  consider  the  question  of  solving  the  projected  Newton  system  (2.8).  It 
is  clearly  best  to  use  an  iterative  method  to  solve  this  linear  system,  so  that 
it  is  not  necessary  to  form  or  store  Z'GZ.   Instead,  a  subroutine  is  required  to 
multiply  Z'GZ  onto  any  vector  y.  This  can  be  accomplished  by  first  multiplying  by 
Z,  then  by  G,  and  then  by  Z~.  Thus  subroutines  are  required  to  multiply  any 
vector  y  by  G,  Z  and  l}   respectively.  The  first  uses: 


-10- 

where  the  terms  corresponding  to  active  residuals^are  omitted  as  explained  in  Sec- 
tion 2.     The  subroutines  which  multiply  a  vector  y  by  Z  and  I^  use  the  list 
structure,  which  represents  Z,  directly.    .These  are  also  required  to  compute  the 
vectors     Z'g  and  Zp^    •     Each  of  the  three  subroutines  requires  0(n)  operations  for 
the  matrix-vector        multiplication.    .:   -  ■  ■ 

Following  [7],  we  do  not  usually  carry  out  the  iterative  solution  of  the  linear 
system  to  full  accuracy.  Let  us  use  the  conjugate  gradient  method  to  solve  the 
system,  denoting  its  iterates  by  w-;  tbe.'relgvsnt- fcnriolas  may  be  found  in  [7], 
with  Wg  =  0.     VJe  terminate     the  i  tc  rati  en, :  setting  py  =  Wj^ ,    if 

IZ  GZw,    +  Z  gll   <.  CGTOLPZ  gH  ,  but  with"  a  maximum' of  k  <_  CGf^AX  imposed     unless 

T 
!Z  gH    <_  CGTHRESH.     This  rule  is   -ather"  arbitrary,  bui  it  allows  avoiding  the 

expense  of  too  many  conjugate  gradient  iterates  wlien  an  accurate  solution  of  (2.8) 
is  not  required,   namely  away  from  the  minimum     on  the  current  manifold.   The  over- 
all work  to  approximately  solve   (2.8)  is'tfius  oT  the  arder  of  n  times  the  number 
of  conjugate  gradient  iterates  used,.. '  "Of'coii"S'e^~,  3't~wbul'd'  be  desirable  to  precon- 
dition    the  conjugate  gradient  it,e^^'tior1,"butri%-j%^r®V  clear  how  to  do  this  since 
Z^GZ  is  not  explicitly  available,     '  '"  -''  -"2"'  ■' '  -    ra   .r.^r^ 

-'-■-£  -i-;-  2sr,  2r/    .dc,-: 

Let  us  now  consider  the  line  se?.rcTi':C2^.S^V^hWe"~  tV  n^  'difficulty  in  using  the 
line  search  of  [14]  unchanged,  sinc^^afV  Bfe  5ieces_sa'ry  ^>,ifonnation  is  available. 
Usually  one  step  of  this  line  sear'ch  ,(|bta;iiis  ;a]  sufficient  decrease  in  F.  This  step 
requires  estimating  the  zeros  of  all  1:'he^  inajqtfive  re^s'i duals  in  the  direction 
along  the  line.  Let  us  say  that  there'. ffCe^CC'tjiactfVe  residuals,  and  let 
{z-},  j  =  l,...,q,  be  the  estimateS'Qf  ^thi'r^eftjs.:  It  ;is  explained  in  [14]  how 
these  values  are  compared  with  an  estimate tQ"^  the  Tnihimum  of  F  along  the  line, 
supposing  that  F  is  smooth.  If  F  dOSs'not_appear'to  be  smooth,  because  of  the 
fact  that  the  latter  estimate  lies  beyond  one'or  more  of  the  (Zj),  it  is  necessary 
to  compute  the  minimum  of  a  piecewise  linear  function  defined  by  the  {z,-}.  This 
computation  amounts  to  finding  a  weighted  median  of  q  points,  and  is  most  easily 
done  by  a  partial  "Heapsort"  of  the  values,  which  requires  0(q  log  q)  operations. 
However,  it  is  possible  to  do  this  computation  in  0(q)  time,  although  this  is  not 
worthwhile  unless  q  is  very  large.  See  [1,15]  for  further  details.  In  any  case, 
this  is  not  a  significant  factor  at  least  in  the  latter  stages  of  minimizing 
(1.3),  since  then  q,  the  number  of  inactive  residuals,  is  much  less  than  m,  the 
number  of  triangles. 

Let  us  turn  our  attention  to  (2.7),  by  far  the  most  difficult  step.  As  already 
mentioned,  the  matrix  A,  when  all  columns  are  included,  is  highly  rank  deficient, 
and  hence  the  Lagrange  vectors  of  (2.4)  are  not  uniquely  defined.  Consequently 
it  is  very  difficult  to  tell  whether  a  given  point  v,  which  is  optimal  on  the 
current  restricting  manifold,  is  in  fact  optimal  in  the  whole  space.  Furthermore, 
even  if  it  v/ere  known  that  v  is  not  optimal,  it  would  be  very  difficult  to  decide 
which  terms  to  delete  from  the  active  set  to  obtain  a  descent  direction.  This  type 
of  difficulty  is  called  degeneracy  in  the  context  of  linear  programming  and 
linearly  constrained  optimization.  The  only  known  methods  guaranteed  to  overcome 
degeneracy  require  time  exponential  in  the  number  of  active  terms,  which  is  out 
of  the  question.  We  therefore  have  the  unsatisfactory  situation  that  we  must 
terminate  the  algorithm  with  a  "solution"  that  we  know  to  be  optimal  on  the  final 
restricting  manifold,  but  for  which  we  are  unable  to  verify  optimality  in  the 
whole  space. 

This  difficulty  seems  to  be  inherent  in  problem  (1.3),  because  of  the  extreme 
built-in  degeneracy.  However,  there  are  a  number  of  positive  remarks  which  can  be 
made.  First,  assuming  that  the  tolerances  of  [14]  are  chosen  small  enough  so 
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that  an  accurate  solution  is  demanded,  note  that  there  is  no  possibility  of  the 
algorithm  terminating  with  too  few  residuals  set  to  zero.  The  only  danger  is 
that  too  many  of  them  may  have  been  set  to  zero.  Second,  the  situation  is 
comparable  to  the  well  known  problem  of  finding  the  global  minimum  of  a  nonconvex 
function.  Several  runs  may  be  made,  using  different  starting  values  or  different 
choices  of  the  parameters,  and  the  lowest  "minimum"  found  may  be  taken  as  the 
best  candidate  for  the  solution.  Third,  and  iiiost  important,  the  numerical 
results  of  the  next  section  indicate  that, eat^least  for  a  certain  initial  guess 
vvO),  the  "solutions"  to  which  the  algorithm  converges  in  practice  are  very   near- 
ly optimal.  Fourth,  we  mention  a  suggestion  of  Dax  [6]  at  the  end  of  the  next 
section.  "5   •  sn:  ,tl. 

We  conclude  this  section  by  desci?ibifig-'anialternative  approach  to  solving  (1.3), 
based  on  ideas  of  [6,8]  and  others:  "This .reqaires  consideration  of  the  hyper- 
bolic approximating  (HAP)  function  of.  [8] ,  giveii.by 

F(v)  -  I     (tr,(v)D2  +  e2)l/2  .  (3.5) 

The  function     F  is  differenti^le^eve^rywlifire  provided  that  e^  >  0.   The  smaller 
the  value  of     e^,  the  better,  .f  ajDproxijTiatJ^s;,  F,     but  a  consequence  of  this,  of   „ 
course,  is  that  the  smaller- ti\€^  v^JupT  o^^^H  ^^  ,t]tfe,' worse  is  the  conditioning  of  F. 
The  alternative  method,  the'n,  simply  requires' mi hiliii zing  F  directly,   using  a 
straightforward  Newton  method.     This  has  the  advantage  that  the  degeneracy  is 
avoided  and  optimality  of  a  fy^i rvt  mijii^i^i net  F^  can  be  verified.     However,   there 
are  four  serious  difficulties  with -thi"?  japp'roach..^    FiVst,  a  minimum  of  F,   not  F, 
is  obtained.     This  may  not  be  ,1-mFkar1^'ri.t»; jsLijic^  0^^.3)   is  derived  from  (1.1)   by  a 
discretization,  but  it  means  ■that  ^^'".sj^ut'i oil i^to  tl.3)   is  not  obtained.   Second, 
if  t^  is  small,  the  extreme  il'L-condi,1^fp;iin^  of  F_makes  the  convergence  of  a 
Newton  method  very  slow.     Thir^J,     eac;b^ffe\^tQn-.ste,l3,  requires  the  solution     of 
a  linear  system  with  much  larger  :dfmensfon^th_cirl'.tliat  of  (2.8),  since  there  are  no 
zero  residuals  to- reduce  the  dimens"ion,  of  .the  .Space.     Especially  in  the  last 
stages  of  the  iteration,  the  dimension 'of  (2.8)   is  much  less  than  n,  and  conse- 
quently (2.8)  can  be  approximately -solved  using  fewer  conjugate  gradient  iterates 
than  a  full   Newton  system  requires.     Fourth,     if     z^     is  small,   the  conditioning 
of  the  full   Newton  system  is  much  worse  than  the  conditioning  of  (2.8),   and  hence 
the  convergence  of  the  "inner"  iteration  can  be  expected  to  be  much  slower  than 
that  for  (2.8).     This  consideration  is  separate  from  the  convergence  of  the 
"outer"  nonlinear  iteration   (the  second  point)     and  the  dimension  of  the  linear 
system  (the  third  point). 

This  alternative  algorithm,  minimizing  F,   is  easily  implemented  in  the  same  pro- 
gram that  was  prepared  for  the  direct  minimization  of  F.      In  particular,   although 
the  alternative  algorithm  does  not  allow  zero  residuals,  the  same  linked  list  data 
structure  is  used  to  impose  the     boundary  conditions  and  the  constraint  (1.3c). 
Thus  the  "full   Newton"  system  is  actually  implemented  in  the  form  (2.8),  but  with 
the  matrix  Z  having  rank  equal   to  the  number_of  variables  minus  the  number  of 
constraints.     This  corresponds  to  a  matrix     A  containing  only  the  columns  which 
describe  the  constraints,   i.e.   the  1x1   blocks  for  the  Dirichlet  conditions,   and 
the  last  column     for  (1.3c). 

NUMERICAL  RESULTS 

The  results  were  obtained  using  Fortran  on  a  VAX  11/780  at  the  Courant  Mathematics 
and  Computing   Laboratory.      Double   precision  arithmetic,   i.e.   aoproximately  16 
decimal   digits  of  accuracy,  was  used  throughout.     We  restricted  the  experiments 
to  Strang's  model    problem,   i.e.    (1.1)  with  f  e  0  and  c  h  1,  whose  solution     is 
shown  in   (1.2).     Because  of  the  symmetry,  we  actually  solved  the  problem  on  the 
lower  left  quarter  of  the  square  shown  in  (1.2).     Thus  the  domain  was  a  unit 
square,  with  N  mesh  points  in  each  direction,  homogeneous  Dirichlet  boundary 
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conditions  on  the  lower  and  left  sides,  and  with  the  points  on  the  other  two  sides 
constrained  to  match  pairwise.  This  is  conveniently  implemented  using  the  linked 
list  data  structure,  with  an  initial  configuration  of  one  fixed  list  for  each 


point  on  the  lower  and  left  sides,  one  variable  list  for 
the  upper  and  right  sides,  and  one  variable  list  for  the 
and  each  interior  point.  "'-'i  c- 


each  pair  of  points  on 
top  right  corner  point 


(0)    •  --J  • 
We  chose  the  initial  vector  v^  '   as  follows,  where 
left  side)  to  N-1  (upper  or  right  side)  r^ '"■-'■!:  :,^ 


i  and  j  range  from  0  (lower  or 


v«"(x,.yj) 


'  0     if  i  =  0  or  j=0  (Dirichlet  boundary  conditions) 

Jl^Oc  "29.'  nruiT'  •:•  .  - 

=  I    N-1    if  i  +  j  setl''-  9rJ  -^D-  r.ic^^   ; 

•B   .-3d ■-'.■;:,.; -jD  ,"  :■ 
i+j-1  otherwise     s-e-;  ,  , 


This  vector  was  then  scaled  to  satisfy  ( 1. 3c )i£."":T^  .-corresponding  piecewise 
linear  function  is  constant,  i.e.  has  zero  nesri<fei£al:S:,^:in  the  shaded  part  of: 


■: ::    1  9\v 


and  is  linear  across  the  remainder  of.  tbe-grid. except  at  the  boundary.  This 
choice  of  v^  '  is  a  compromise  between  v^/=  constant  in  the  interior,  which  has 
too  many  zero  residuals,  and  v=  linear,  not  constant,  in  the  interior,  which 
having  no  zero  residuals  causes  Z'GZ  to  be  singular  initially. 


The  results  are  summarized  in  Tables  I  and  II. 
solving  (1.3)  when  high  accuracy  is  demanded. 
were  set  as  follows:  Em^u  =  10"^°,  c.^t  =  10' 


N  =  33' 


'LINE 


=  10 
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-MCH 


=  0.9.  Of 


'ACT 
these, 


Table  I  gives  the  results  of 

In  this  case  the  parameters  of 

znr   =  10"°  (except  10"^  for 

which 


[14] 


'PG 


the  most  important  is 

is 


z^g. 


to  be  made, 


The 


specifies  how  small  the  norm  of  the  projected  gradient,  ..  <_  3..,  .^  ^^  ^^   ...^^w.  ..... 

parameter  tj^Qj   determines  when  an  attempt  is  made  in  (2.6)  to  make  "almost"  active 
residuals  "exactly"  active.  This  parameter  may  be  reduced  dynamically  by  the 
program;  see  [14]  for  details.  The  conjugate  gradient  iteration  parameters  were 
set  as  follows:  CGTOL  =  10"',  CGTHRESH  =  10"^,  and  CGMAX  as  shown  in  the  table. 
The  other  columns  in  the  table  specify  N  (the  number  of  mesh  points  in  each  direc- 
tion), the  column  rank  of  Z  at  the  computed  solution  (i.e.  the  dimension  of  the 
final  restricting  manifold),  the  final  value  of  F,  ITER  (the  number  of  major 
iterates  required),  and  CGITER  (the  total  number  of  conjugate  gradient  iterates). 
The  most  interesting  thing  to  note  in  Table  I  is  the  consistency  of  the  final 
values  of  F  determined  by  the  various  choices  of  CGMAX.  The  dependence  of  CGITER, 
the  best  measure  of  the  total  amount  of  work,  on  CGMAX,  is  also  interesting 
although  somewhat  inconclusive.. 

The  minimizing  solution  v  is  shown  for  the  case  N  =  13  in  Table  III.  Note  that 
there  are  two  regions  where  v  is  constant,  the  large  one  in  the  top  right  and  the 
small  one  in  the  bottom  left.  The  band  of  varying  values  between  the  two  regions 
is,  of  course,  the  result  of  the  discretization.  One  can  see  that  the  circular 
arc  r  of  (1.2)  is  reflected  in  this  band  of  varying  values  since,  for  example, 
v(x2,y2)=  0.209  is  greater  than  v(xi,y3)  =  v(x3,yi)  =  0.149. 
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In  order  to  show  how  well  the  solution  of  (1.3)  approximates  the  solution  of  (1.1) 
as  N  is  increased,  we  display  in  Figure  1  graphs  of  the  solution  v  plotted  along 
the  diagonal  of  the  square.  The  results  from  Table  I  are  shown  for  N  =  9,  17 
and  33,  The  step  function  shown  is  the  analytic  solution  of  (1.1).  The  distance 
from  the  mesh  point  (0,0)  to  the  discontinuous  jump  in  the  solution,  divided  by 
the  length  of  the  diagonal ,  is  o(l  -  1  /  'II)    =  0.155,  where  -p  is  given  following 
(1.2).  The  analytic  solution  has  a  value  u-.=-l '/rtl  -  p  +  7Tp'^/4)  =  1.054.  in  the 
constant  positive  region,  with  a  corresponding  minimal  F  value  of  2-2p+Trp/2 
=  1.772. 

.-! -:'■-  j '  Z  -\ 
Table  II  summarizes  results  when  a  much  less  accurate  solution  of  (1.3)  is 
required.  Here  we  show  results  both  for  the  direct  minimization  of  F  and  for  the 
minimization  of  the  HAP  functiori  F,  described  at'the  end  of- Section  3.  In  the 
latter  case,  we  set  eu  =  10"  -(h^/Z)  ,  where  h  =  1/(N-1),  since  the  residual 
sizes  are  proportional  to  the  areas  of  the  triangles.  If  a  much  larger  value  of 
EH  is  used,  the  function  F  ':app.r-oxiiTnates  F  Very  poorly;  if  a  much  smaller  value 
is  used,  the  computation  time  ijecSDntes 'exjcessive.  We  set  epg ,  the  tolerance  on 
the  projected  gradient,  to  10"  .  Again*  a  "luch  smaller  value  of  cpg  leads  to 
excessive  computation  time,  since  F  js  ill -conditioned.  We  therefore  also  used 
epQ  =  10"'  for  the  runs  which  minimi 2£_.E-<i4rertTy,^.l^e  also  set  e^j-y  =  10"^  for 
these  runs,  except  c/^qj  =  10"  for  N  ■  33T""tyiis  js  actually  relevant  only  for 
N  <_  9,  since  for  larger  values  of  itL-tb€-resT3u'al  sizes  are  smaller.  The  other 
parameters  had  the  same  values  as  for^Tatle-  f.  Note  that  it  is  F,  not  P,  which 
is  shown  in  the  column  for  the  HAP 'resu.lis,  as  well  as  for  the  direct  results. 
The  interesting  observation  to  make  TrQirf Table  II  is  that  even  when  an  inaccurate 
solution  is  required,  the  method  which  directly  minimizes  F  has  substantial 
advantages  over  the  HAP  method.  Iir  order  to  show  the  difference  in  the  accuracy 

of  the  solutions  from  Tables  I  and  IIj the  solution  obtained  from  minimizing 

the  HAP  function  when  N  =  17  is  plotted' in  Figure  IT 

Dax  [6]  has  made  an  interesting  suggestion' regarding  the  use  of  the  HAP  function 
in  the  context  of^multi facility  location  problems.  He  suggests  switching  to  the 
minimization  of  F  only  after  reaching  the  minimum  of  F  on  the  current  manifold, 
as  a  way  of  potentially  escaping  from  a  nonoptimal  point.  If  a  lower  value  of 
F  is  obtained  in  the  process,  a  switch  can  be  made  back  to  the  minimization  of  F. 
This  seems  a  useful  idea,  especially  if  v(0)  is  poorly  chosen,  e.g.  v^^'  =  con- 
stant in  the  interior.  In  this  case  thg  direct  algorithm  would  temiinate 
immediately,  but  a  switch  to  minimizing  F  produces  a  lower  value  of  F  without 
difficulty.  However,  this  idea  does  not  produce  any  benefits  for  the  runs  given 
in  Table  I,  apparently  because  the  direct  minimization  produced  nqv^j   nearly 
optimal  results. 

CONCLUDING  REMARKS 

In  this  paper  we  have  confined  our  attention  to  the  description  of  a  method  which 
obtains  accurate  solutions  to  (1.3),  a  discretization  of  (1.1).  Some  alternatives 
to  (1.1)  which  we  have  not  considered  include  allowing  discontinuous  elements 
(Johnson  [11])  and  solving  the  primal  problem  for  which  (1.1)  is  the  dual  (Strang 
[16]).  Other  relevant  papers  include  [5,10,13,17,18].  We  conclude  with  the 
remark  that  if  a  method  such  as  the  one  described  in  this  paper  is  to  be  used  to 
obtain  very  accurate  solutions  to  (1.1)  or  related  problems,  some  sort  of  mesh 
refinement  near  the  curve  of  jump  discontinuity  would  probably  be  desirable. 
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TABLE  I 
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TABtg  II  or:: 

Results  of  Minimizing  Both  F  and  the  HAP 


.,Ujg:=10i^ 
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TABLE  III 
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